
import Analysis 

#[ ..., [r,m,Num_bosons,sigma,Num_stars],...]
# Args = [
#     [0.5,1.0,0,1,10000],
#     [0.5,1.0,10000,1,10000],
#     [1,0.5,20000,1,10000],
#     [5,0.1,100000,1,10000],
#     [10,0.05,200000,1,10000],
#     [50,0.01,1000000,1,10000],
#     [0.5,1.0,10000,1,0]
# ]

Args = [
    [1,1,0,0.001,1000],
    [1,1,0,0.0002,5000],
    [1,1,0,0.0001,10000],
    [1,1,0,0.00002,50000],
    [1,1,0,0.00001,100000]
]

for args in Args:
    print("----------New Analysis--------")
    print(
        f"r = {args[0]}",
        f"mu = {args[1]}",
        f"Num_bosons = {args[2]}",
        f"sigma = {args[3]}",
        f"Num_stars = {args[4]}"
    )
    Analysis.analysis(*args)

/home/boris/Documents/Research/FDM_n_Bodies/1D_Codes/Non-Dim/Analysis
----------New Analysis--------
r = 1 mu = 1 Num_bosons = 0 sigma = 0.001 Num_stars = 1000

/home/boris/Documents/Research/FDM_n_Bodies/OneD/WaveNonDim.py:129: RuntimeWarning: invalid value encountered in true_divide
  F_s = F_s/Norm_const

v_rms = 0.5923277051957598
z_rms = 0.15871675883738892
K_avg = 0.5*m*v_rms^2 = 0.1754260551712375 (m=1)
=> 2*K_avg = 0.350852110342475
W_avg = 158.7167588373889
------------------
K_tot = 0.17542605517123733
K_avg = 0.00017542605517123733
W_tot = -0.14596981745731766
W_avg = -0.00014596981745731765

/home/boris/Documents/Research/FDM_n_Bodies/1D_Codes/Non-Dim/Analysis/Analysis.py:266: RuntimeWarning: invalid value encountered in true_divide
  v_rms_array = bins/bins_counts

/home/boris/Documents/Research/FDM_n_Bodies/1D_Codes/Non-Dim/Analysis/Analysis.py:442: RuntimeWarning: divide by zero encountered in log
  ax[1].plot(np.log(z_right),np.log(rho_whole))

Check
Check
Check

----------New Analysis--------
r = 1 mu = 1 Num_bosons = 0 sigma = 0.0002 Num_stars = 5000

v_rms = 0.5851737606787921
z_rms = 0.1541549084183391
K_avg = 0.5*m*v_rms^2 = 0.17121416509348014 (m=1)
=> 2*K_avg = 0.3424283301869603
W_avg = 770.7745420916955
------------------
K_tot = 0.1712141650934794
K_avg = 3.424283301869588e-05
W_tot = -0.14558193639607853
W_avg = -2.9116387279215706e-05

Check
Check
Check

----------New Analysis--------
r = 1 mu = 1 Num_bosons = 0 sigma = 0.0001 Num_stars = 10000


import matplotlib.pyplot as plt
v_rms_s = [0.5923277051957598,0.5851737606787921,0.581649317774426,0.5731593297761615,0.5840606962146646]
z_rms_s  = [0.15871675883738892,0.1541549084183391,0.15731951911022882,0.15992388602013577,0.1588461256252678]
Num_p_s = [1000,5000,10000,50000,100000]

fig,ax=plt.subplots(1,2,figsize= (15,5))
ax[0].plot(Num_p_s,z_rms_s,'o--')
ax[0].set_title("$z_{rms}$ vs Number of Particles")
#ax[0].set_ylim(0,0.5)

ax[1].plot(Num_p_s,v_rms_s,'o--')
ax[1].set_title("$v_{rms}$ vs Number of Particles")

plt.show()


import numpy as np
import matplotlib.pyplot as plt

My_Package_PATH = "/home/boris/Documents/Research/Coding"
import sys
sys.path.insert(1, My_Package_PATH)
import OneD.NBody as NB

z = np.linspace(-1,1)

x = 0
v = 1
star = NB.star(0,1,x,v)

dt = 0.1
t = 0
i = 0
while t < 2:
    plt.plot(star.x,0,'ro')
    plt.xlim(-1,1)
    plt.show()

    star.x -= v*dt
    star.reposition(2)
    t += dt


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import LogNorm, Normalize
import os
import subprocess
import cv2 
from PIL import Image 
import scipy.optimize as opt

#Import My Library
My_Package_PATH = "/home/boris/Documents/Research/Coding"
import sys
sys.path.insert(1, My_Package_PATH)
import OneD.WaveNonDim as ND
import OneD.NBody as NB
import OneD.GlobalFuncs as GF

#Set up Directory for saving files/images/videos
# Will not rename this again
dirExtension = "1D_Codes/Non-Dim/Analysis"
Directory = os.getcwd()#+"/"+dirExtension #os.curdir() #"/home/boris/Documents/Research/Coding/1D codes/Non-Dim"
print(Directory)

r,m,Num_bosons,sigma,Num_stars = [0.5,1.0,0,1,10000]

mu = m #M_scale = 1

L = 2
N = 10**3
z = np.linspace(-L/2,L/2,N)
dz = z[1]-z[0]

folder = "ParticlesOnly_Snapshots"
stars_x = np.loadtxt(folder+"/"+f"StarsOnly_Pos.csv", dtype = float, delimiter=",")
stars_v = np.loadtxt(folder+"/"+f"StarsOnly_Vel.csv", dtype = float, delimiter=",")
Energies = np.loadtxt(folder+"/"+"Energies.csv", dtype = float,delimiter = ",")
#chi = np.loadtxt(folder+"/"+f"Chi.csv", dtype = complex, delimiter=",")
chi = np.zeros_like(z)
centroids = np.loadtxt(folder+"/"+"Centroids.csv",dtype = float, delimiter=',')

stars = [NB.star(i,sigma,stars_x[i],stars_v[i]) for i in range(len(stars_x))]

grid_counts = NB.grid_count(stars,L,z)

rho = (grid_counts/dz)*sigma 


i = 0
max_bool = False
while max_bool == False:
    for j in range(len(rho)):
        if rho[j] > rho[i]: #if you come across an index j that points to a larger value..
            #then set i equal to j
            i = j 
            #break
        else:
            max_index = i
            max_bool = True
max_rho = rho[max_index]

#Other method to accumulate left and right sides:
for star in stars:
    star.x = star.x - z[max_index] #shift
    star.reposition(L) #reposition

grid_counts = NB.grid_count(stars,L,z)
rho_part = (grid_counts/dz)*sigma 
#Add the density from the FDM
rho_FDM = mu*np.absolute(chi)**2 
rho = rho_FDM + rho_part

centroid_z = 0
for j in range(len(grid_counts)):
    centroid_z += z[j]*grid_counts[j]
centroid_z = centroid_z / Num_stars

stars_x = [star.x for star in stars]

std = np.std(stars_v)
mean_x = np.mean(stars_x)


R = 0
while True:
    R += dz
    mass_enclosed = 0
    star_collection = []
    for star in stars:
        if np.abs(star.x-mean_x) <= R:
            mass_enclosed += 1
            star_collection.append(star)
    print(R,mass_enclosed)
    if mass_enclosed >= 0.5*Num_stars:
        break

print(R)
plt.figure()
plt.scatter(stars_x,stars_v,s=1)
xx = np.linspace(-R,R,100)
plt.plot(xx,np.sqrt(R-xx**2))
plt.plot(xx,-np.sqrt(R-xx**2))
plt.scatter([star.x for star in star_collection],[star.v for star in star_collection],s=1)
plt.show()

Sigma = std**2 / R
print(Sigma)

/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis

---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb Cell 5' in <cell line: 37>()
     <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000004?line=34'>35</a> stars_x = np.loadtxt(folder+"/"+f"StarsOnly_Pos.csv", dtype = float, delimiter=",")
     <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000004?line=35'>36</a> stars_v = np.loadtxt(folder+"/"+f"StarsOnly_Vel.csv", dtype = float, delimiter=",")
---> <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000004?line=36'>37</a> Energies = np.loadtxt(folder+"/"+"Energies.csv", dtype = float,delimiter = ",")
     <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000004?line=37'>38</a> #chi = np.loadtxt(folder+"/"+f"Chi.csv", dtype = complex, delimiter=",")
     <a href='vscode-notebook-cell:/home/boris/Documents/Research/Coding/1D_Codes/Non-Dim/Analysis/CurveFitting.ipynb#ch0000004?line=38'>39</a> chi = np.zeros_like(z)

File /usr/lib/python3/dist-packages/numpy/lib/npyio.py:1067, in loadtxt(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin, encoding, max_rows, like)
   1065     fname = os_fspath(fname)
   1066 if _is_string_like(fname):
-> 1067     fh = np.lib._datasource.open(fname, 'rt', encoding=encoding)
   1068     fencoding = getattr(fh, 'encoding', 'latin1')
   1069     fh = iter(fh)

File /usr/lib/python3/dist-packages/numpy/lib/_datasource.py:193, in open(path, mode, destpath, encoding, newline)
    156 """
    157 Open `path` with `mode` and return the file object.
    158 
   (...)
    189 
    190 """
    192 ds = DataSource(destpath)
--> 193 return ds.open(path, mode, encoding=encoding, newline=newline)

File /usr/lib/python3/dist-packages/numpy/lib/_datasource.py:533, in DataSource.open(self, path, mode, encoding, newline)
    530     return _file_openers[ext](found, mode=mode,
    531                               encoding=encoding, newline=newline)
    532 else:
--> 533     raise IOError("%s not found." % path)

OSError: ParticlesOnly_Snapshots/Energies.csv not found.


G = 6.67E-11
print(R)
print("--------------------")
print("")

Sigma = std**2 / (np.pi* R**(3/2))
print(Sigma)

print(10000/R)
print(std**2)
print(10000/(np.pi*R**2))

print(std**2 * R)

print("--------------------")
print("")
new_std = np.std([star.v for star in star_collection])
Sigma = new_std**2 / (np.pi* R**(3/2))
print(f"Sigma = {Sigma}" )

print(10000/R)
print(new_std**2)
print(10000/(np.pi*R**2))

print(new_std**2 * R)

0.048048048048046965
--------------------

100502.33917739333
208125.0000000047
3325.3676414862753
1378791.600352967
159.77742421555317
--------------------

Sigma = 61907.311466298816
208125.0000000047
2048.3560084912815
1378791.600352967
98.41950791549479


v_rms = np.sqrt(np.mean([star.v**2 for star in stars]))
print(v_rms)

57.66620191650019


v_mean = np.mean([star.v for star in stars])
std = np.sqrt(np.sum([(star.v - v_mean)**2 for star in stars])/(len(stars)-1))
print(std)

Sigma = std**2 / (np.pi* R**(3/2))
print(Sigma)

print(10000/R)
print(std**2)
print(10000/(np.pi*R**2))

print(std**2 * R)

57.66888425752163
100512.39041643498
208125.0000000047
3325.7002115074265
1378791.600352967
159.79340355590878


phi_part = GF.fourier_potentialV2(rho_part,L)
phi_part = phi_part - np.mean(phi_part)
print(np.mean(phi_part))

phi_part = phi_part - np.max(phi_part)

# Compute Chandrasekhar's potential energy tensor:
a_part = NB.acceleration(phi_part,L)
W = 0
for i in range(len(z)):
    dW = rho_part[i]*z[i]*a_part[i]
    W += dW
print(W)

a_part = NB.acceleration(phi_part,L)
W = 0
for i in range(len(z)):
    dW = -0.5*rho_part[i]*phi_part[i]
    W += dW
print(W)

# Compute only for the stars that exist:
a_part = NB.acceleration(phi_part,L)
W = 0
for star in stars:
    g = NB.g(star,a_part,dz)

    dW = - star.x*g
    W += dW / Num_stars
print(W)

# phi_part = GF.fourier_potentialV2(rho_part,L)
# a_part = NB.acceleration(phi_part,L)
# W = 0
# for i in range(len(z)):
#     dW = - dz*a_part[i]**2 / (8*np.pi)
#     W += dW
# print(W)
#W = np.sum(phi_part)
#print(W)

# Compute only for the stars that exist:
W = 0
for star in stars:
    #g = NB.g(star,a_part,dz)
    i = int(star.x//dz)
    rem = star.x % dz 

    if i != len(phi_part)-1:
        value = phi_part[i] + rem*(phi_part[i+1]-phi_part[i])/dz
    elif i == len(phi_part)-1:
        # then i+1 <=> 0
        value = phi_part[i] + rem*(phi_part[0]-phi_part[i])/dz
    
    phi_star = value
    dW = phi_star
    W += dW
print(W)

0.0
-16785655529.943058
55851483242.48459
-1500.1026550328993
-7621006.55368937


# Compute total KE of stars:
K = 0
for star in stars:
    dK = 0.5*sigma*star.v**2
    K += dK
print(K)
#average KE:
print(K/Num_stars)

# #Compute Total Potential
# W = 0
# for star in stars:
#     #g = NB.g(star,a_part,dz)
#     i = int(star.x//dz)
#     rem = star.x % dz 

#     if i != len(phi_part)-1:
#         value = phi_part[i] + rem*(phi_part[i+1]-phi_part[i])/dz
#     elif i == len(phi_part)-1:
#         # then i+1 <=> 0
#         value = phi_part[i] + rem*(phi_part[0]-phi_part[i])/dz
    
#     phi_star = value
#     dW = phi_star
#     W += dW
# print(W)
# #average W:
# print(W/Num_stars)

# Compute only for the stars that exist:
a_part = NB.acceleration(phi_part,L)
W = 0
for star in stars:
    g = NB.g(star,a_part,dz)

    dW = - sigma*star.x*g
    W += dW
print(W)
print(W/Num_stars)

16626954.217372933
1662.6954217372934
-15001026.55032902
-1500.102655032902


v_rms = np.sqrt(np.mean([star.v**2 for star in stars]))
z_rms = np.sqrt(np.mean([star.x**2 for star in stars]))
print(f"v_rms = {v_rms}")
print(z_rms)
#v_rms = np.sqrt(np.sum([star.v**2 for star in stars])/Num_stars)

K = 0.5 * v_rms**2
print(f"K_avg = 0.5*m*v_rms^2 = {K} (m=1)")
print(F"=> 2*K_avg = {2*K}")

print(z_rms*Num_stars)

print("-----------------------")




R_syst = Num_stars / v_rms**2
print(R_syst)



rho_0 = np.mean(rho_part)
print(4*rho_0*z_rms)

print(v_rms**2 / (2*np.pi*z_rms))

print(16*np.pi*rho_0**2*z_rms**3 / Num_stars)

v_rms = 57.66620191650019
0.1567842541769929
K_avg = 0.5*m*v_rms^2 = 1662.6954217372852 (m=1)
=> 2*K_avg = 3325.3908434745704
1567.842541769929
-----------------------
3.0071653140030277
3132.549398456389
3375.673107160588
483.3349205244898


plt.plot(z,phi_part)
plt.plot(z,-Num_stars/np.abs(z))
plt.ylim(5*np.min(phi_part),-np.min(phi_part))

(-126523.1647523825, 25304.632950476498)


phi_part = phi_part - (np.max(phi_part)-np.max(-Num_stars/np.abs(z)))

plt.plot(z,phi_part)
plt.plot(z,-Num_stars/np.abs(z))
plt.ylim(5*np.min(phi_part),-np.min(phi_part))
plt.show()

# Compute total KE of stars:
K = 0
for star in stars:
    dK = 0.5*star.v**2
    K += dK
print(K)
#average KE:
print(K/Num_stars)

#Compute Total Potential
W = 0
for star in stars:
    #g = NB.g(star,a_part,dz)
    i = int(star.x//dz)
    rem = star.x % dz 

    if i != len(phi_part)-1:
        value = phi_part[i] + rem*(phi_part[i+1]-phi_part[i])/dz
    elif i == len(phi_part)-1:
        # then i+1 <=> 0
        value = phi_part[i] + rem*(phi_part[0]-phi_part[i])/dz
    
    phi_star = value
    dW = phi_star
    W += dW
print(W)
#average W:
print(W/Num_stars)

16626954.217372933
1662.6954217372934
-107621006.55368945
-10762.100655368946


def f(z,*p):
    u_0 = p[0]
    z_0 = p[1]
    return u_0 / np.cosh(0.5*z/z_0)**2

guess = [rho_0,z_0]
popt,pcov = opt.curve_fit(f,z,grid_counts,p0 = guess)
plt.plot(z,grid_counts)
plt.plot(z,f(z,*popt))
plt.show()

guess = [rho_0,z_0]
popt,pcov = opt.curve_fit(f,z,phi_part,p0 = guess)
plt.plot(z,phi_part)
plt.plot(z,f(z,*popt))
plt.show()

def g(z,*p):
    return p[0]*np.exp(-z**2 / p[1])

guess = [-rho_0,z_0]
popt,pcov = opt.curve_fit(g,z,phi_part,p0 = guess)
plt.plot(z,phi_part)
plt.plot(z,g(z,*popt))
plt.show()

Compute Total KE and Total Potential Energy of Stars¶

Calculate $v_{rms}$ and $R_{syst}$¶